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Abstract The Ornstein-Uhlenbeck process may be used to generate a noise signal with a finite corre- 
lation time. If a one-dimensional stochastic process is driven by such a noise source, it may be analysed 
by solving a Fokker-Planck equation in two dimensions. In the case of motion in the vicinity of an 
attractive fixed point, it is shown how the solution of this equation can be developed as a power series. 
The coefficients are determined exactly by using algebraic properties of a system of annihilation and 
creation operators. 
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1 Introduction 

Many physical processes are modelled by adding noise to a dynamical system. This paper discusses 
a one-dimensional example where noise is added to a system with an attracting fixed point. With a 
suitable change of variables the fixed point is positioned at the origin and the system may be expressed 
in the form 

x = -x- eg{x) + f{t) (1) 
where f{t) is a random noise with statistics 

</(*)> = 0, (f(t 1 )f(h)) = C(t 1 -t 2 ) (2) 

(throughout this paper (X) denotes the expectation value of any random variable X). In {T]) the non- 
linearity of the stochastic process is represented by the function g(x), which is assumed to satisfy 
g{0) — <?'(0) = 0, with a multiplier e which will be used as a perturbation parameter. In the case where 
f{t) is a white noise signal, with correlation function C{At) = 2D 8 {At), the probability density for 
the solution of ([1]) satisfies a Fokker-Planck equation: 

dP d d 2 P 

_ = _ [(l + e9W) P 1+0 _. (3) 

In one dimension it is easy to obtain exact steady state solutions of this equation. When e = 0, it 
is also possible to determine the propagator and use this to determine correlation functions exactly 
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(this case is known as the Ornstein-Uhlenbeck process [T], described in [HE]). The exact results which 
are available when e = can be used as a basis for a perturbation expansion in e (an example of the 
application of this approach is described in |4J). When the noise is not delta-correlated however, it is 
much more difficult to analyse equation Jl} . This paper considers a particular case where the noise in 
([1]) is not delta-correlated, but where it is nevertheless possible to analyse the statistics of the solutions 
by means of a Fokker-Planck equation. This is possible if the noise is itself generated by a dynamical 
process which is driven by a white noise signal. 
The following model will be analysed: 

x = —x — eg(x) + y 

y = -uy + (i) (4) 

where is a white-noise signal with statistics 

(£(*)> = 0, (t{t 1 )£(fr)) = 28(t l -t 2 ). (5) 

The variable y(t) is not influenced by x(t). It is a coloured noise signal which is generated by an 
Ornstein-Uhlenbeck process, which has an exponential correlation function [IJEIE] 

(y(*lM*2)) =wexp{-u\ti -t 2 \) ■ (6) 

This approach to modelling systems with coloured noise generated by an Ornstein-Uhlenbeck process 
was previously considered by Fox et al [5], who pointed out that the numerical simulation of systems 
with coloured noise is facilitated by using an Ornstein-Uhlenbeck process to generate the noise, and 
by Risken [3J, who noted that it allows process with coloured noise to be modelled using a Markovian 
process. The process represented by equations (j4j, §5§ and ((6|) approaches the process described by ((3]) 
in the limit as u> — -> oo, with the diffusion constant D — 1 (a scaling of x and t can always be applied 
so as to set D = 1). 

In many applications the precise form of the correlation function of the noise is not known, and the 
analysis of this model is sufficient to understand the effect of the noise having a finite correlation time. 
In cases where the precise form of the correlation function of the noise is significant, the approach 
which is developed here can be extended to more general correlation functions by generalising the 
dynamical process which is used to smooth the white-noise signal. 

The joint probability density for x and y satisfies the Fokker-Planck equation 

BP f) f) f> 2 P 

^ = j^-y^m +U ^P) + ^ w (7) 

which will be written as 

^ = [P + eT l ]P (8) 

(a 'hat' over a symbol will be used to denote a differential operator). Our objective is to obtain the 
solution of the steady-state Fokker-Planck equation in the form of a power series 

oo 

P(x,y) = J2^Pj(^y) ■ (9) 

3=0 

Inserting this into (JSj) gives the recursion relation Pj+i(x, y) = —!FQpxPj(x,y). It is not immediately 
clear how the inverse P^ 1 can be computed or whether this equation gives a meaningful expression 
for Pj+i. The operator P^ 1 will be defined by analysing the spectrum of J-q. It is clear that there is 
an eigenfunction of jFq with eigenvalue equal to zero, which is the steady-state solution of |(7J| when 
e = 0. This suggests that the inverse T^ 1 is ill-defined. However, it will be shown that the terms 
in the series are in fact well-defined. Moreover, the coefficients can be determined exactly by using 
the algebraic properties of a system of raising and lowering operators, which are defined in section^ 
before developing the perturbation theory in section [3] The use of annihilation and creation operators 
to treat Fokker-Planck equations is discussed in the book by Risken [3], but the usual approach is 
not applicable to the problem which is treated here, and a different method is required. This point is 
considered in section 0J 
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A motivation for developing this approach was to study the Lyapunov exponents for inertial parti- 
cles suspended in random fluid flows. This problem can be transformed into determining the expectation 
value of a stochastic variable such as (Q}, where f(t) represents the velocity gradient of the fluid at the 
position occupied by the particle 013] • The technique which is developed here will be applied to the 
calculation of the Lyapunov exponent in a companion paper j7. . In this work, however, the method 
will be developed in a general setting, which will surely find applications in other areas. 



2 Unperturbed eigenfunctions and spectrum 

The steady state solution of the Fokker-Planck equation ((7|) for the e = problem is a Gaussian 
function 

1 r(l + uj) ((I + uj)x 2 + y 2 ~ 2xy) 



P {x,y) = exp 



2w 2 



(10) 



Writing P — /Pq, where P is an eigenfunction of the Fokker-Planck equation satisfying TqP = XP, it 
is found that / satisfies 



(x - y)d x f + [2(1 + u)x - (2 + uj)y]d y f + ufd$f = F f = Xf 



(11) 



The operator F defined in (|f 1 1) maps polynomials to polynomials, and the generalised eigenfunctions 
of Fo satisfying Fof(x,y) = Xf(x,y) must be polynomials in x and y. Clearly fo(x,y) = 1 is an 
eigenfunction with A = and by inspection the linear functions fi(x,y) = (I + uS)x — y and f u (x,y) = 
2x — y have eigenvalues A = — 1 and A = — lo respectively. This observation motivates the definition 
of 'lowering operators', B\ and B^, which satisfy B\f^ = and = 0. Such operators can be 

constructed as linear combinations of derivatives, d x — d/dx and d y = d/dy. By inspection, it is 
possible to construct a set of raising and lowering operators which satisfy the commutation relations: 





= -A l 


Fq, 


Fo, Bi 


= B X 


FqjBu 



(12) 



u)B u 



(where [A, B] = AB — BA). These relations imply that the eigenvalues are A nm = — (n + Lorn). The 
required operators are 



Au. 
B, 



0.r 



2x-y-^-r 
9 X + (1 + (jj)d y 



u>{u)-l) a 

uj+i u v 



A 1 = (I + d)x -y- ^jd x 



(13) 



B x =d x + 2d, 



v • 



The operators A\ and A u are termed raising operators because they increase the quantum number n 
and m respectively. The operators B\ and B u are termed lowering operators. 

Rather than working with the polynomials which are eigenfunctions of fo it is more convenient to 
work with eigenfunctions of the operator fio defined by 0, (jHJ), satisfying fo<j>nm = -(fi+wm)^, m . The 
operators in (|I3p can be transformed into raising and lowering operators for generating eigenfunctions 
of To directly. These are 



-d x + (uj - l)d y 



01 = ^1 [d x + 2d y ] +x+^- iy K = t^ + (1 + ^Kl - zhv 
These satisfy the commutation relations 



(14) 



&i,0i 



-1 



-1 



(15) 
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and 



-0L\ 



3~o, Pi — (3\ 



(16) 



The inverse relations to (fTl)) are: 
X = -&i 



x = $i + + ^ [wdi - o,„] y = (1 - + [2wai - (w + l)da 
Expressing .Fq in terms of the raising and lowering operators yields 

Fb = - fwd w &, + axflt 



(17) 



(18) 



The notation will be used as a shorthand for a function ip(x, y) (this is an adaptation of the Dirac 
notation of quantum mechanics). The eigenfunctions <p n ,m{x,y) of jFq are denoted by vectors \4> n ,m) 
and their eigenvalues are given by 



Fq \4>n, m ) = -(n + Lom) \<j> n>m ) 



(19) 



Because Tq is not self-adjoint, these eigenfunctions are not orthogonal. Neither are they assumed to be 
normalised. The commutation relations (| 1 5[) . (|16p for the raising and lowering operators are consistent 
with the following relations defining how the eigenfunctions may be generated by successive application 
of the raising operators: 



ai\4>n,m) = \4>n+l,m) &u>\<f>n,m) = \4>n,m+l) 

Pl\4>n,rn) = Tl \4> n ~l,rn) $w\<l>n,m) = ™ \<t>n,m-l) ■ 



(20) 



The first two relations enable all other eigenstates to be generated by successive application of the 
raising operators starting from the steady state |</>o,o)- It will be assumed that |0o,o) is normalised as 
a probability density, but the other eigenstates need not be normalised. 



3 Iteration of perturbation series expansion 

Consider the iteration of the perturbation series expansion ([9]) starting from \Pq) = |0o.o) ■ Assuming 
that each term is expanded in terms of the eigenfunctions: 

oo oo 

fi) = EE C nllW. (21) 
n— m— 

The following considers how the perturbation series may be evaluated for the case where g{x) is ex- 
pressed as a polynomial in x. The case g(x) — x 2 is considered explicitly, but more general polynomials 
are treated in exactly the same way. 

Successive terms in the series expansion are given by 

\Pj+i) = -f^Fi \Pj) = ^ai&lPj) ■ (22) 

The definition of the inverse Fq 1 appears problematic, because one of the eigenvalues of J-q is equal to 
zero. Note however that the operator F^ 1 acts on a state which is multiplied by the raising operator 
&i, so there is no division by zero. More explicitly, if a function Q{x,y) is expressed in the form 



oo oo 

Qn,m $n,m ( x i y) 

n — rn—0 



(23) 
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then, using (j2"0|) . 



.X. V _^ 



•^0 ^llQ) = - E E 9n,m — T 1 - |0n+l,m) (24) 



n— m— 

^2 



which is well-defined. The operator x may be expressed in terms of raising and lowering operators via 
(fTT)) . so that x 2 \Pj) can be expressed as a linear combination of the form (|23[) with known coefficients. 

It is, therefore, possible to iterate to determine the coefficients c™ in (|2ip . starting from Cnm = &no8 m o- 
The successive corrections \Pj) are generated from (f2"2"| . using (fT7]), (f^Uj) and (|2"4"|) . 

It is often moments of the distribution which are required, rather than the probability density 
itself. These can be obtained from the expansion (12ip if the moments of the eigenfunctions \4> n ,m) are 
known. First, note that the raising operators are constructed from partial derivatives. It follows (using 
integration by parts) that, except for the case n = m = 0, 

/CO 
dy 4> n , m (x,y) = S n<0 S m . (25) 
-oo 

To determine other moments of the functions \4> n ,m,), such as (x N ), express x and in terms of raising 
and lowering operators using (|17[) . so that x N (f> nym (x,y) may be expressed as a linear combination of 
eigenstates: 

oo oo 

x -EE 

n'—0 m'—O 



)=E E W(N 7 n,m,n',m')\<f> n ,, m ,) (26) 



where the coefficients W(N,n,m,n' , to') are readily obtained from the commutation relations of the 
raising and lowering operators. Upon integration over x and y the only contribution comes from the 
term where n' = ml = 0, so that 

/oo />oo 
dx / dy x N <f> n , m {x, y) = W{N, n, to, 0, 0) . (27) 
-oo J — CO 

As an example, consider the evaluation of (x)\ 

2 ^ 

X\4>nm) = ™|</>n-l,m) +«l|(/> nm _i) + — 5 |0 n +l,m) + ~, ?\4>n.rn+l) (28) 

CtT — 1 1 — LO z 

so that the only cases where W(l, n, to, 0, 0) ^ are W(l, 1, 0, 0, 0) = 1 and W(l, 0, 1, 0, 0) = 1. This 
means that if the probability density P{x, y) is expressed in the form (|23p with coefficients p n , m then 
(x) = pifi +po.i- From (f2"4"]l it can be seen that the application of Pq a.\ to \Q) yields a state for which 
the coefficient <?o,i is equal to zero. It follows that the series expansion of (x) is 

oo 

<*> = E eic S- (29) 

3=0 

The coefficients can be determined using an algebraic manipulation package. They are equal to zero 
for even values of j. The first three non-zero coefficients are 



(3) _ w 2 (10cj 2 + 27w + 15) 

~ (cc + l) 3 (l + 2^) 
(5) _ 10 (72 oj 6 + 540 oj 5 + 1630 w 4 + 2493 w 3 + 2010 oj 2 + 807 w + 126) w 3 

~ (1 + 2w) 2 (u; + 1) 5 (1 + 3c 1 j) (2 + w) 



(30) 



The complexity of the coefficients increases rapidly as j increases. In the limit as lu — > oo, where the 
coefficients approach those for the case of a white noise signal, the coefficients are already known. As 
lu — > oo the non-vanishing coefficients approach 1,5,60, 1105, . . ., which is in agreement with setting 
r = in the results contained in 0]. 
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4 Concluding remark 

The technique described here may find applications in a wide variety of systems which can be modelled 
by stochastic differential equations. A companion paper will describe the application of the method to 
the calculation of Lyapunov exponents for particles moving in a random fluid flow. 

One technical point shown be remarked upon. There is a standard approach to linear Fokker-Planck 
equations in which the Fokker-Planck operator To is transformed to a Hermitian operator TCq, writing 

Ho = f-^ot (31) 

where the transformation T is a 'gauge transformation', exp['P(x, y)]. By a suitable choice of a phase 
function <P(x, y) which is a quadratic form in x and y the non- Hermitian contributions to Tio can usually 
be eliminated. This approach is discussed in [3j. The Hermitian form of the operator is convenient for 
subsequent calculations, because its eigenfunctions are orthogonal. Moreover, because Tio is a quadratic 
form in position and momentum operators Xi and pi = —id Xi respectively, Ti.o can be transformed into 
a harmonic oscillator Hamiltonian, which allows standard harmonic oscillator raising and lowering 
operators to be used directly (this approach was used in [4]). This method is, however, not applicable 
to the operator To defined by ([7]) and l[B]). because the non-Hermitian components of H.q cannot be 
eliminated by any choice of the coefficients of the quadratic form (p(x,y). 

There are two possible ways to avoid this difficulty. One is to explore whether the required results 
can be obtained using the set of non-orthogonal eigenfunctions of To generated by the operators 6l\ 
and a u . This is the approach which has been adopted here. It has been shown that series expansions 
of (a;^) may be obtained using an expansion in terms of eigenfunctions which are neither orthogonal 
nor normalised. 

In some other contexts it may still be desirable to transform To to a Hermitian operator, so that 
orthonormal bases can be used. It is possible to make a more general transformation of the form ([3"Tj) . 
which is closely related to the definition of the Weyl representation of metaplectic operators in quantum 
mechanics, as discusses in [5]. To this end it is useful to define an operator 

/oo 
dX exp(-X 2 /2a)exp{-Xd x ) (32) 
-oo 

and a similar operator A4 y (a), in which d x is replaced by d y . To understand the significance of 
this operator, consider its action upon xf(x). In the following manipulations exp(— Xd x ) is identi- 
fied with its Taylor series and hence interpret it as a translation operator which shifts a function by X: 
cxp(—Xd x )f(x) — f(x — X). Assuming that f{x) is a normalisable and sufficiently smooth function, 



M x xf(x)= / dX exp(-X 2 /2a)(x-X)f(x- X) 

J — oo 

/oo 
dX exp(-X 2 /2a)[xexp(-Xd x )f - Xexp(-Xd x )f] 
-oo 

= xM x f + a 



f°° d 

J dX — [exp(~X 2 /2a)} exp(-Xd x )f 

f°° d 
= xM x f -a J dX exp(-X 2 /2a)— exp(-Xd x )f 

= xM x f + ad x M x f . (33) 

Because M x (a) commutes with d x , (|33|) implies two equivalent rules for commuting x and M. x (a): 

M x {a)x = (x + ad x )M x {a) , xM x {a) = M x {a){x - ad x ) . (34) 

Note that M. x {a) also commutes with d y , and y. The action of operators M x {a) and Xi y ((3) enables 
J-q to be converted into a partial differential operator which is of second order in both x and y. It is 
found that the application of the operator Ai x (a) was not required to convert To to self-adjoint form, 
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however it is necessary to introduce a multiplication by a scalar function. A transformation operator 
of the form 

f = M y (l3)exp(-$) (35) 

can be used in (j3l"j) . where $(x, y) is a quadratic form <&{x, y) = \{Ax 2 + By 2 + 2Cxy). The parameters 
A, B, C, and (3 can be chosen so that Hq is self-adjoint. 

Acknowledgement. I am grateful to Michael Morgan for a careful reading of the manuscript. 
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